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Abstract 

The requirements for accurate numerical simulation are increasing constantly. Modern high precision physics ex- 
periments now exceed the achievable numerical accuracy of standard commercial and scientific simulation tools. One 
example are optical resonators for which changes in the optical length are now commonly measured to 10^^^ precision 
|[T|. The achievable measurement accuracy for resonators and cavities is directly influenced by changes in the distances 
between the optical components. If deformations in the range of 10~^^ occur, those effects cannot be modeled and 
analysed any more with standard methods based on double precision data types. New experimental approaches point 
out that the achievable experimental accuracies may improve down to the level of 10^^*^ in the near future 0. For the 
development and improvement of high precision resonators and the analysis of experimental data, new methods have 
to be developed which enable the needed level of simulation accuracy. Therefore we plan the development of new 
high precision algorithms for the simulation and modeling of thermo-mechanical effects with an achievable accuracy 
of 10~^°. In this paper we analyse a test case and identify the problems on the way to this goal. 



Motivation 

Optical high-finesse resonators are widely used as a 
frequency reference for the stabilization of lasers e.g. 
in optical atomic clocks, now aiming for the 10~^^ 
level of precision [3], or in direct tests of special and 
general relativity. For example, ID and ||4l tested 
the isotropy of the speed of light by comparing the 
frequencies of two lasers referenced to crossed res- 
onators with a relative precision of 10-^5 at 1 s inte- 
gration time. 

Typically, an electronic feedback loop is applied to 
actively stabilize the laser frequency to the optical 
resonator. The laser frequency then results directly 
from the geometrical resonator dimensions as given 
by vl = mc = 2L, where L is the resonator length 
and m is the longitudinal mode number. The mea- 
surement performance thus depends on the length 
stability of the implemented cavities. Since Av/v = 
AL/L, high precision frequency stabilization at the 
10^^^ level and below, requires a relative cavity 
length stability of the same magnitude. The domi- 
nant effects which cause expansions or contractions 
are gravity (in particular if an experiment is assem- 
bled on earth and conducted in space) and thermal 
expansion of the components. Therefore high preci- 
sion modeling of thermo-mechanical effects is a ba- 
sic requirement for a safe design and verification of 
optical resonators. Only if the models and simula- 
tions offer the same precision which can be achieved 
in experiments, the optimization, design, and eval- 
uation of optical resonators can be performed in a 



sound way. 

To date, the best cavity stabilized lasers have achieved 
a relative frequency stability on the order of few parts 
in 10^^ |I3. As pointed out by Numata et al. 161 
|71, this precision is limited by fundamental thermal 
noise, mostly affecting the highly reflective mirror 
coatings. Consequently, low-noise, single-crystalline 
silicon mirrors are currently being developed |2| and 
an improvement of optical cavity stability down to 
the 10^^ level might become feasible already within 
the next years. These cavities will ultimately help to 
build optical clocks of improved accuracy or to carry 
out fundamental tests at the 10^^^ level and below. 
This development however needs to be paralleled by 
similar advancements in high precision modeling 
tools. 

In order to react to the grown achievable accuracies 
of optical resonators and the need for improved mod- 
eling methods, we intend to develop a new method 
for the modeling of thermo-mechanical effects in 3 
dimensions with an achievable numerical accuracy 
of 10~^'^. In this paper, we outline the difficulties to 
overcome by analysing a cavity test case. 

Test case configuration 

In the following we will present the cavity test case 
which will be used to demonstrate the problems for 
modeling of thermo-mechanical effects with a nu- 
merical accuracy of lO^^'^. The analysis is confined 
to a single mirror within the cavity and performed 
using two different approaches. 
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Fig. 1 : Design of cavity test case configuration. 



Cavity test case 

An optical cavity is typically axially symmetric and 
composed of two highly reflective concave lenses 
which are bonded to a spacer. The dimensions of the 
test cavity analyzed here are adapted from Webster 
et al. iSI, for the design see Figure [T] The mir- 
ror substrates, the mirror coatings, and the spacer 
are assumed to be made of ultra-low expansion glass 
(ULE) provided by Corning with the relevant mate- 
rial parameters listed in Tab. [T] Although ULE is 
known for anisotropic thermal expansion, we assume 
as an approximation that all materials are homoge- 
nous and isotropic and that the space around the cav- 
ity is a vacuum at room temperature. Further, we take 
only thermo-elastic effects into account and neglect 
other error sources like vibration or material losses. 



Material parameter 


value 


Thermal conductivity [W/mK] 


1.31 


Specific heat [J/kgK] 


767 


Mean coefficient of thermal ex- 


± 30 X IQ-^ 


pansion [1/K] 




Elastic modulus [Pa] 


67.6 X 10^ 


Poisson's ratio 


0.17 


Mass density [kg/m'^] 


2.21 X 10^ 



Tab. 1 : Material parameters of ULE Coming code 7972. 



space and cannot be deformed. As the mirrors are 
bonded to the spacer, a deformation of the mirrors at 
the contact areas is also not allowed. However, heat 
conduction at the contact area is taken into account. 
ULE is designed to have a minimum thermal expan- 
sion coefficient at some temperature Tule near room 
temperature, whose exact value has to be determined 
experimentally for each cavity. Even if Tule is ex- 
actly known, in practice it is not possible to achieve 
a stable temperature of Tule over the whole cav- 
ity, and a small temperature gradient is inevitable. 
However, by careful thermal control a temperature 
stability of a few mK can be obtained if operating 
near Tule and even about 100/LiK for temperatures 
around 300K [9|. Here we assume the mirrors sub- 
strates and coatings to have a stable temperature cor- 
responding to a coefficient of thermal expansion of 
30 X 10~^K"^. For more rigorous treatments than 
this test case the thermal expansion coefficient should 
be varied in time to reflect the small disturbances of 
the temperature. 

Also, we assume that this test cavity has a finesse of 
F = 10^ and a transmission coefficient of the incou- 
pling mirror of T = 10^^. For an input power of 



Br 



Modeled effects, boundaries, and loads 

Regarding the setup of the test cavity, the relevant 
thermal physical effects are heat conduction and ra- 
diation as well as the thermal expansion of the ma- 
terial. Due to the assumed room temperature, heat 
radiation may be neglected as a first simplification. 
The spacer may radiate to the whole surrounding space 
and, thus, we assume it as a heat sink with a con- 
stant temperature of 300K as a second simplifica- 
tion. Third, we assume that the spacer is fixed in 



bOfiW this results with the gain factor of 

T2 



— T 



10^ 



(1) 



in an effective power of about P = 5W stored in the 
cavity. The losses A due to absorption corresponding 



to the chosen values for F 



T+A 



and T can be 



assumed to A ^ 10^^ resulting in a dissipated power 
of about 5fiW per mirror. The laser is assumed here 
to have a diameter of 200/im and to hit the mirror 
perfectly at its center. 

Once an equilibrium state is reached, the tempera- 
ture is constant in time but may vary over the cavity. 
However, due to the changed temperature the cavity 
will deform, what in turn changes the boundary con- 
ditions for the heat conduction. Therefore, the final 
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static state has to be determined iteratively. For this 
test case, we neglect the influence of the changes of 
boundary conditions and perform one static analysis, 
only. 

FE models and results 

The test case described above was modeled using the 
Finite Element (FE) approach. Two independ pro- 
grams were used for this task: the academic version 
of the commercial ANSYS software and the software 
package MaiProgs developed at the If AM. 

ANSYS model 

As a first approach, one half of the axially symmetric 
cavity was modeled in three dimensions. However, 
due to the limited node number of 256000 nodes of 
the academic version of ANSYS, this turned out to 
be nonsufficient for reaching a high modeling ac- 
curacy. Therefore, the axial symmetry of the prob- 
lem was used to reduce the model to two dimen- 
sions. For the computation of heat conduction and 
mechanical deformation, the eight nodes PLANE77 
and PLANE 183 element types were used. On the 
spacer a fixed temperature and no displacement are 
prescribed and, therefore, we modeled it in ANSYS 
with a uniform mesh with a large maximal element 
size of 1mm. On the lens we started with the same 
element size and reduced it step by step. However, 
the best results could be obtained by using a mesh 
adapted to the problem with maximal element sizes 
of approximatly 0.002mm in the regions of interest 
and a total node number of about 247000. The max- 
imum temperatures for the various degrees of free- 
dom are shown in figure [2l After solving the ther- 
mal problem the resulting temperature distribution 
was applied as a body load for the following struc- 
tural analysis. From a preliminary analysis the dis- 
placements were expected in the range of 10~^^mm. 
Therefore, as the structural analysis is a linear prob- 
lem, we decided to rescale the thermal expansion co- 
efficient for a better numerical performance. The re- 
sulting displacements at front and back of the cavity 
are shown in figures [3] and 01 

MaiProgs model 

At the IfAM computations were made for linear el- 
ements on tetrahedral meshes with up to seven mil- 
lion nodes. For the computations the self developed 
FE/BE-software MaiProgs |[TOI was used. In Figures 
|5] and [6] the axial displacements on the cross section 
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Fig. 2: Max. temp vs. DOF using ANSYS 
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Fig. 3: Nodal solution plot using ANSYS 

of the cavity at the front and the back is plotted using 
the data of the computations with MaiProgs. 
To point out the sensitivity of the mathematical model 
on the accuracy of the solution, computations on two 
different meshes were performed at the IfAM. In the 
first case the usual h- version on a quasiuniform mesh 
was computed. The geometry is modeled via two 
cylinders where the ring of the outer cylinder fits 
with the ring of the cavity that is bonded to the spacer. 
The tetrahedral elements have nearly the same shape 
and size. Since the laser is striking the cavitiy in 
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Fig. 4: Nodal solution plot using ANSYS 
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Fig. 5: Nodal solution plot using MaiProgs 
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Fig. 6: Nodal solution plot using MaiProgs 

its midpoint with a radius of 0.1mm it is obvious 
that the temperature will have a large gradient there. 
This is also emphasized by the axial displacement 
in Figure [5] where the displacement seems to have 
a singularity in the midpoint of the cavity. Using 
this information the IfAM modeled a second mesh 
consisting of three cylinders where the inner cylin- 
der has the radius of the laser such that the coars- 
est mesh has 400 more points in the neighborhood 
of the midpoint than the first version. The tetrahe- 
dral elements do not have the same size any longer, 
since the elements near the midpoint of the cavity are 
considerably smaller than the other elements. From 
this point of view the second mesh can be regarded 
as kind of an adaptive mesh. This leads to a bet- 
ter approximation of the boundary loads and hence 
to a better approximation of the solution as can be 
seen in Figure |7] Here the red dashed line is the ex- 
trapolated maximal temperature inside the cavity. On 
both meshes the IfAM computed four steps of the h- 
version ITTl [T2l starting with a coarsest mesh with 
2400 points and 2800 points, respectively. Figure |7] 
shows that the computations on the second mesh lead 
to much better approximations of the temperature. It 
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Fig. 7: Max. temp vs. element size 

seems reasonable to assume that an adaptive compu- 
tation using some kind of a posteriori error estimator 
(residual, hierarchical, etc.) llT3l IT4l would exceed- 
ingly improve the accuracy of the solution. 

Comparison of results 

The two different softwares and approaches described 
in this section lead to results comparing well to each 
other. However, a closer inspection of the two results 
reveal significant differences: the maximal tempera- 
tures for the two results differ by 

At = 1.904041846501059 x 10~°^ K . 

Given this quite large discrepancy it is almost supris- 
ing that the maximal displacements for both approaches 
differ by only 

Ad = 1.447725122807002 x 10"^^ mm. 

These differences are caused on the one hand by the 
discretization errors connected to the different meshes, 
and on the other hand by the limited accuracy of the 
numerical procedures involved in the solution pro- 
cess. Whereas the first error source can be mini- 
mized by using finer meshes, the second one is an 
inert problem of floating point arithmetic. As a rule 
of thumb, one can say that for a result with n sig- 
nificant accurate digits an accuracy of 2n significant 
digits for each addition and multiplication is needed. 
As the exponent is treated seperatly, for linear prob- 
lems with known order of magnitude for displace- 
ment or temperature the accurate digits can be con- 
siderably improved by a suitable reseating of the ma- 
terial properties. However, for nonlinear problems 
this approach will in general not work. As a con- 
sequence, the only way to minimize the numerical 
errors is to use quad (with 128 bits) or even arbitrary 
precision arithmetic instead of the usual double pre- 
cision (with 64 bits). Since there is no appropriate 
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Fig. 8: Thermal energy error norm plot 

hardware available, these computations must be per- 
formed with a software emulation. This slows down 
the computation process significantly. 

Accuracy Evaluation 
ANSYS 

Within the ANSYS software a posteriori error esti- 
mator can be used to evaluate the quality of the so- 
lution. As the exact solution is not known in gen- 
eral and, in particular, not for our test case, this error 
estimator compares the (discountinous) heat flux or 
stresses, respectively, in each node as computed by 
ANSYS to an averaged (continuous) function, which 
interpolates these values. In each element Ei the en- 
ergy norm of the difference of the discontinuous 
values and the interpolated function is a measure for 
the error in this element. For the whole model a per- 
centaged relative energy norm error can be given by 
the sum of all Cj's divided by Cj and the thermal 
dissipation energy. This is a kind of Zienkewics/Zhu 
error estimator |[T5l . which converges to the true er- 
ror. The percentage error for various degree of free- 
doms are plotted in Figures [8] and |9]for both the ther- 
mal and structural computations. 



Fig. 9: Structural energy error norm plot 
relative errors in energy norm 



Mai Progs 

Since the exact solution of the presented test case is 
not known in advance it is impossible to compute 
the exact error of the computed approximations in 
the energy norm. Therefore Figure in [TO] extrapo- 
lated norms for the temperature and the displacement 
are compared with the energy norms of the approxi- 
mated solutions. The error does not reflect the accu- 
racy of the solution but can be considered as a mea- 
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sure to compare different approximations. The aim 
is to derive efficient and reliable a posteriori error 
estimators for the test case in order to have a better 
prediction of the true error. Of course, these estima- 
tors can be used within adaptive algorithms leading 
to a faster convergence 1.13. 141 . 



Conclusion 

We analysed a cavity test case for thermal induced 
changes in the cavity length, which directly influ- 
ences the stability of the laser frequency used to in- 
terrogate the cavity. We used two different and inde- 
pendent softwares for a finite element analysis of the 
thermo-elastic deformation of the cavity lens due to 
the dissipated power of the laser. 
Despite the numerous simplifications assumed for this 
test case the results of the two independently com- 
puted values for the maximal displacement differ by 
an order of magnitude of 10^^^, which is consid- 
erably lower than the corresponding achievable fre- 
quency stabilities of ultrastable cavities. The individ- 
ual analysis have two major error sources: the dis- 
cretization error of the mesh and the numerical error 
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of the solution process. The first error can be shown 
to converge to zero for an infinitly fine mesh. Succes- 
sively finer meshes can be produced with any com- 
mercial FE software and adequate computing power. 
Therefore, this error source can be minimized with- 
out principle problems. The second error source is 
inert to floating point arithmetic and can only be min- 
imized if data types with more bits are used. For 
the present hardware, these data types have in gen- 
eral to be software emulated, which will slow down 
the computation process considerably. Also, to the 
knowledge of the authors, a publicly available FE- 
solver using quad or arbitrary precision arithmetics 
does not exist. 

The test case presented here is processed taking into 
account only thermo-elastic effects. For the design 
of high precision cavities, this is of course not suffi- 
cient. However, the principle problems for the FE- 
modelling will remain the same. Therefore, cavi- 
ties with length stabilities on the 10^^'' level or even 
higher can only be simulated if the problem of nu- 
merical errors can be solved. We plan to solve this 
problem with the development of a stand-alone FE- 
solver using arbitrary precision arithmetics, which is 
based on the solver already implemented in MaiProgs, 
together with an interface to ANSYS or other pre- 
pocessors. 
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